namespace NBody.KD
{
    using System;
    using NBody;
        
    public class SplitNode : Node
    {
        private int cut_dim;
        private double cut_val;
        private Node left;
        private Node right;
        
        public SplitNode(int d, double p, Node left, Node right, double[,] bounds)
        {
            this.cut_dim = d;
            this.cut_val = p;
            this.left = left;
            this.right = right;
            
            this.n = left.NumParts + right.NumParts;
            this.mass = left.Mass + right.Mass;
            this.cm = new double[3];
            for (int i = 0; i < 3; i++)
                this.cm[i] = left.Mass / this.Mass * left.CentreOfMass[i] +
                             right.Mass / this.Mass * right.CentreOfMass[i];
                
            ComputeRmax(bounds);
        }
            
        public override void FindNearest(double rd, ref SearchInfo info)
        {
            int cd = cut_dim;
            double old_off = info.Off[cd];
            double new_off = info.Position[cd] - cut_val;
            if (new_off < 0)
            {
                left.FindNearest(rd, ref info);
                rd += -old_off*old_off + new_off*new_off;
                if (rd < info.PQ.TopPriority)
                {
                    info.Off[cd] = new_off;
                    right.FindNearest(rd, ref info);
                    info.Off[cd] = old_off;
                }
            }
            else
            {
                right.FindNearest(rd, ref info);
                rd += -old_off*old_off + new_off*new_off;
                if (rd < info.PQ.TopPriority)
                {
                    info.Off[cd] = new_off;
                    left.FindNearest(rd, ref info);
                    info.Off[cd] = old_off;
                }
            }
        }
        
        public override void CalcPotential(IDParticle p, double theta, double eps, ref double phi)
        {
            double eps2 = eps * eps;
            double x = p.X - cm[0];
            double y = p.Y - cm[1];
            double z = p.Z - cm[2];
            double r = Math.Sqrt(x * x + y * y + z * z + eps2);
            double r_open = 2.0 / Math.Sqrt(3.0) * r_max / theta;
            
            // if particle is within r_open, open cell and kepp computin'
            if (r < r_open)
            {
                left.CalcPotential(p, theta, eps, ref phi);
                right.CalcPotential(p, theta, eps, ref phi);
            }
            else // just add in cell's c/m
            {
                phi -= mass / r;
            }
        }
    }
}
